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ACCURACY OF FINITE ELEMENT APPROXIMATIONS 
TO STRUCTURAL PROBLEMS 

By Joseph E. Walz, Robert E. Fulton, Nancy Jane Cyrus, 
and Richard T. Eppink 
Langley Research Center 

SUMMARY 

This paper reports on a theoretical investigation of the convergence properties of 
several finite element approximations in current use and assesses the magnitude of the 
principal errors resulting from their use for certain classes of structural problems. 

The method is based on classical order of error analyses commonly used to evaluate 
finite difference methods. Through the use of Taylor series the differential or partial 
differential equations which represent the convergence and principal error characteris- 
tics of the finite element equations are found. These resulting equations are then com- 
pared with known equations governing the continuum and the error terms are evaluated 
for selected problems. Finite elements for bar, beam, plane stress, and plate bending 
problems are studied as well as the use of straight or curved elements to approximate 
curved beams. The results of the study provide basic information on the effect of inter- 
element compatibility, unequal size elements, discrepancies in triangular element 
approximations, flat element approximations to curved structures, and the number of 
elements required for a desired degree of accuracy. 

INTRODUCTION 

♦ 

Finite element methods have been used for many years with good success in the 
analyses of complex structures and many aerospace structures are designed on the basis 
of these analyses. Although the use of these methods is widespread, not enough is known 
of the theoretical accuracy and convergence properties of finite element models when 
used to represent structures. Accuracy studies are usually based on numerical solutions 
to specific problems for comparison with known results. Convergence studies are car- 
ried out by investigating the convergence of the numerical results as the number of ele- 
ments is increased. (See ref. 1, for example.) Such methods, while valuable in providing 
a cursory assessment of the adequacy of various model approximations, are heavily 
dependent on the numerical data and the problem studied and may obscure the true 


character of the approximation. More basic information is needed on the accuracy and 
convergence properties of various finite elements for structural approximations. 

The purpose of the present paper is to report on an analytical investigation of the 
accuracy of several finite element stiffness method approximations in current use and to 
assess the magnitude of the errors resulting from the use of these approximations for 
certain classes of structural problems. The present paper includes results which were 
summarized in references 2 and 3 and in addition provides a comprehensive explanation 
of the background details of the error analysis procedure which led to these results. 

The method used is based on classical order of error analyses commonly used to evaluate 
the discretization errors of finite difference methods. Through the use of Taylor series 
the ordinary or partial differential equations are found from which the convergence and 
error characteristics of the finite element equations can be determined. These resulting 
equations are compared with the known continuum equations governing the structure. 

The discretization error in the finite element approximation is evaluated for a limited 
class of deflection and vibration problems to provide simple formulas for the size of an 
element required to obtain a certain degree of accuracy. (See refs. 4 and 5.) Finite 
elements for bar, beam, plane stress, and plate bending problems are studied as well as 
the use of straight and curved elements to approximate a curved beam. 

SYMBOLS 


A 

a, b 
B 


c 

D 

E 

F 

Fx, Fy 
h 


cross-sectional area of one -dimensional element 

length of rectangular panel in x- and y-directions, respectively 

Et 


extensional stiffness of plate, 
constant 

bending stiffness of plate, 
Young’s modulus 
finite element nodal force 


1 - ii* 


Et* 


12(l 




nodal force in x- and y-directions, respectively, of plane stress element 


reference length of finite element 


2 



I 

u 

M 

L 


Mx, My 
m, n 
m 


N 

O(h) 

P 

Po 

P x > Py 

q 

q 

<*o 

R 

t 

u, v, w 


u 0 


moment of inertia of cross section 

i,jth grid point 

element stiffness matrix 

length of one -dimensional structure 

bending moments in x- and y-directions, respectively (see fig. 7) 

harmonic wave numbers for sinusoidally distributed loads 

mass per unit length for one -dimensional element; mass per unit area 
for two-dimensional element 

number of elements per harmonic wave number, or - — 

hm hm 

omitted discretization error term of order h 

tangential loading on bar and arch structures 

amplitude of sinusoidal tangential or inplane loading 

distributed inplane loading in x- and y-directions, respectively 

transverse loading 

magnitude of uniform pressure 

amplitude of sinusoidally distributed transverse loading 
radius of arch 
thickness of plate 

displacements in x-, y-, and z -directions, respectively 
amplitude of sinusoidally varying displacement 



x, y, z rectangular coordinates 

a constant indicating ratio of element dimensions 

vector of nodal displacements 

e discretization error in displacement or frequency parameter 

e j, eg discretization errors in displacement and frequency parameters, respec- 
tively, for plate bending examples 

e , e y , e w discretization errors in displacements u, v, and w, respectively 
6, cp nodal rotation variables for bending element about y- and x-axes, respectively 

I u. Poisson’s ratio 

p radius of gyration of cross section 

ct x , cry, T X y normal stresses in x- and y-directions and shear stress, respectively 
oj circular frequency 

Subscripts: 
ex exact 

i identification for ith grid point 

Abbreviation: 

ACM Adini, Clough, and Melosh plate bending model 

Prime or Roman numeral with a symbol denotes differentiation with respect to x. 

Subscript following a comma denotes differentiation with respect to the indicated 
variable. 


4 



ERROR ANALYSIS PROCEDURE 


Two main sources of error result from the use of finite element methods to solve 
structural problems. These may be conveniently separated into round-off error and 
discretization error. Round-off error is that error associated with the accuracy with 
which numbers are manipulated in a computer and is not considered in the present paper. 
Discretization error is that error associated with using discrete variables to represent 
a problem where the state variables are continuous. This error occurs irrespective of 
the accuracy of numerical calculations and occurs in structural problems when finite 
elements are used to approximate a continuous structure. Discretization errors may be 
of two kinds, as follows: (1) Errors which are a function of the size of the element and 
which vanish as the element size vanishes; and (2) errors which do not vanish when the 
element size vanishes. The relative merit of elements and/or patterns which lead to the 
first kind of discretization errors depends on the relative rate at which the error vanishes 
with element size. Elements and/or patterns which lead to the second kind of discretiza- 
tion errors are unsatisfactory approximations and should be recognized and avoided. 

The present paper deals with an assessment of both kinds of discretization errors in 
finite element approximations. 

The method used in the study is to obtain the typical finite element equations which 
express force equilibrium at a reference node point in terms of displacement variables. 
These finite element equations (which are a class of difference equations) are then 
expanded in Taylor series about the nodal point to obtain the differential equations equiv- 
alent to the finite element equations at that node. The resulting differential equations 
are compared with the governing equations for the continuum approximated. A simple 
bar element approximation is treated in detail as an example to characterize the method 
and to define the terms to be used in the study. 

Example Discretization Error Analysis 

The force-displacement relations of a typical one-dimensional structural element 
having ends i - 1 and i are 



where <{F^> and <0Q> are the vectors of the nodal forces and displacements at the ith 

node and [k] is the element stiffness matrix. Consider a bar of constant cross- 
sectional area A subjected to a distributed axial load p(x) and approximated by 
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finite elements where x denotes distance along the bar (fig. 1). For a bar finite ele- 
ment is the nodal extensional force, is the corresponding displacement uj, 

and 




( 2 ) 


where h is the length of the element and E is Young's modulus. 


In finite element methods distributed loads are replaced by concentrated loads at 
the node points. There are several possible ways in which distributed loads may be 
converted to concentrated loads. In this paper a very simple lumping procedure for 
loads is used for all finite element models. This procedure for the bar gives the con- 
centrated load as the value of the distributed load at a grid point multiplied by one-half 
of the total length of the two adjoining elements. Treating the distributed load in this 
manner and utilizing equations (1) and (2) leads to the following equation for equilibrium 
at an ith point located between two segments of length h and ah (fig. 1) 


EA 

h 




p.h(l + a) 

H ■ _J— 2 


(3) 


Equation (3) is a typical finite element equilibrium equation for the indicated 
approximation. 

The behavior of the system of equation (3) is investigated as the number of equa- 
tions approaches infinity and the size of the element vanishes. This is done by assuming 
the displacements ui to be analytical functions and by examining equation (3) in the 
limit as h approaches zero with the aid of Taylor series expansion of displacements at 
points i - 1 and i + 1 about the ith point. This procedure results in 


u" - |(1 - a)u'" 


+ 


hVl + a^\ iv 

12 Vl + a) 


+ 



(4) 


where the primes denote differentiation with respect to x and where the subscript i 
has been omitted from Uj. (Omission of such subscripts is done consistently throughout 
this paper.) Equation (4) is the differential equation equivalent to the finite element 
equation at the ith node. This is a differential equation of infinite order and a detailed 
discussion of the solution of such equations is given in reference 6. The treatment of 
this equation herein is somewhat heuristic but satisfactory for the present results. The 
terms in equation (4) which are not multiplied by powers of h comprise exactly the 
governing differential equation for the continuous bar and the remaining terms are the 
discretization error in the differential equation that result from use of the finite element 
approximations. Thus equation (3) reduces to the correct continuum equation as h 
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approaches zero. The principal error in the differential equation resulting from the 
finite element approximation is the set of terms in the discretization error containing 
the lowest power of h. The power of h in the principal error denotes the order of 
the discretization error of the finite element equation and the rate of convergence as 
h vanishes. Thus the finite element equation leading to equation (4) has a discretiza- 
tion error of order h. Note that if the segments are equal (a = 1) , the discretization 
error is of order h 2 . Similarly, if the finite element equations do not converge to the 
governing differential equations, the discretization error would be of order h® or 1. 


Harmonic Loading and Vibration Example 

The error analysis procedure described in the previous section gives only the rate 
of convergence of the finite element approximation. The magnitude of the discretization 
error in the bar finite element approximation is now evaluated for the special case of a 
bar supported at each end, subjected to a sinusoidally distributed static loading 

P = P 0 sin 

and approximated by equal length elements. In equation (5) p Q is the amplitude of the 
loading, m is the number of half waves, and L is the length of the bar. For this case 
equation (4) becomes 

1.2 . P~ 

( 6 ) 


( 5 ) 


u 


- . h" „iv 


4- — — u 
12 


° S in2E=0 


EA 


A solution to this equation provides explicitly the error in the finite element approxima- 
tion as a function of element size and harmonic wave length. If the loading is regarded 
as a continuous function, as was done in references 4 and 5, the exact solution to equa- 
tion (6) is 


u = u 


o 


sin 


mux 

L 


(7) 


where u D is the amplitude of the displacement. Substitution of equation (7) into equa- 
tion (6) results in 


u 


u = 


ex 


1 u ex(l + e ) 


1 - e 

where the exact solution to the bar for the sine loading is 


u ex 


_ p o ( L \ 2 . mux 
EA Vmu/ s L 


the principal error in deflection due to the finite element approximation is 

,2 


e = — 


12N 2 


( 8 ) 

( 9 ) 


( 10 ) 
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( 11 ) 


and the number of elements per harmonic half wave used to approximate the bar is 

N = L / m 
h 

Error results for various values of N are as follows: 


N 

percent 

3 

10 

4 

5 

9 

1 


Thus approximately three elements per deflection half wave are needed to keep the error 
in finite element deflection calculations at a node to within 10 percent and nine elements 
per half wave to limit the error to 1 percent. 


A similar approach is used to determine the error in natural frequency of the bar 
example when approximated by equal length elements. For vibration behavior where a 
simple lumping process analogous to that used for the distributed loading is used to 
obtain a diagonal mass matrix, the counterpart of equation (6) is 


u"+^u iv + 


12 


mw 2 

EA 


u = 0 


(12) 


where m is the mass per unit length and a> is the circular frequency. The vibration 
mode shape is the same as equation (7) and the eigenvalues of equation (12) are 


where 


co 2 = J ex (l - e) 



(13) 

(14) 


and e is the discretization error due to the finite element approximation given by equa- 
tion (10). Note that the error in the square of frequency is the same as that resulting 
from static deflection calculations except of opposite sign. Thus the frequency calcula- 
tions converge from below and the deflection calculations converge from above. 

For completeness, table I summarizes the principal error terms for the general 
bar finite element approximations and table II summarizes the evaluation of these errors 
for the harmonic response examples. 
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RESULTS AND DISCUSSION 


The method described in the previous section was used to investigate the conver- 
gence of beam elements, straight and curved element approximations to an arch, rec- 
tangular and triangular plane stress elements, and plate bending elements. The results 
of the investigation are summarized in this section together with a general discussion. 

A summary of results of the error studies is given in tables I and II. 


One -Dimensional Bending Element 

For beam bending problems the well-known finite element model is based on the 
nodal variables 



(15) 


where Wj is the displacement and the rotation at the ith node. (See fig. 2.) The 
governing differential equation for a beam is 

w iv = i (16) 


where q is the lateral load, I is the moment of inertia of the beam cross section, and 
the constraint relationship between rotation and displacement is 

e - w’ = o (17) 

A satisfactory finite element approximation should converge to equations (16) and (17). 

As shown in appendix A (see also table 1(a)), the beam element does converge with a 
principal error term of order h 4 . 

The principal error was evaluated for a finite element approximation to a simply 
supported beam of length L subjected to sinusoidally distributed lateral load q where 

q = q Q sin 2M* (18) 

and where q Q is a constant. The principal error in the lateral deflection w of the 
beam resulting from the finite element approximation is (see table II) 


e = 


720N 4 


(19) 


Equation (19) also gives the magnitude of the error in frequency determination resulting 
from the finite element approximation. (See table II.) 
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Approximation of Curved Structures 


Straight elements are often used to approximate curved structures such as curved 
beams, arches, and shells. To gain some insight into the influence of curvature, an arch 
of radius R was approximated by conventional straight beam elements which also had 
extensional capability. The arch loading is a normal pressure q and a tangential dis- 
tributed loading p. (See fig. 3.) At a typical ith node three variables are required to 
define element behavior. For this problem it is convenient to take these variables as 



(20) 


where u and w are the tangential and radial displacements, respectively, and 9 is 
the rotation. Three finite element equations result from force and moment equilibrium. 
As the element size vanishes, the moment equilibrium equation at the ith node converges 
to the correct constraint equation between the rotation 6 , radial displacement w, and 
tangential displacement u 

I) - o (21) 

and the tangential and normal equilibrium equations converge, respectively, to 


(w' 




( 22 ) 


♦£)-»-" (23 > 

where the principal error terms for equations (21) to (23) are proportional to h 2 . In 
equations (22) and (23) the rotation 6 has been eliminated by using equation (21) in a 
manner similar to that done previously for the other bending problems. Equations (22) 
and (23) result when the cylindrical shell equations given in reference 7 are specialized 
to the case of an arch. Thus, the straight elements provide a convergent approximation 
to a first approximation arch theory with an error of order h4 

A study was also made of the use of curved elements to approximate a curved struc- 
ture. The stiffness matrix for a curved element was derived based on the strain energy 
for a ring given in reference 8 and the details are given in appendix B. The displace- 
ments were approximated by assuming that arch tangential and normal displacements 
were linear and cubic, respectively, over the curved element length. The resulting finite 
element equations were investigated and the element pattern was found to converge to 
equation (21) plus the following equations with errors of order h 2 (see table 1(a)) 


10 





-p-0 


(24) 


E t iv+2 M +EA (-* + §)--° < 25) 

The arch equations (24) and (25) are consistent with the ring strain energy given in 
reference 8 and are also the ring equations which result from the specialization of the 
cylindrical shell theory given in reference 9. Thus the curved elements also provide a 
convergent approximation to the curved structure with an error of order h^. The two 
sets of arch equations (22), (23) and (24), (25) are related in that the curvature terms 
differ by terms composed of the extensional strain divided by the arch radius. According 
to the Koiter criterion for thin shells (ref. 10) such modifications are admissible alterna- 
tives for a first approximation theory of shells, and it seems appropriate to use this 
criterion for arches. 

An assessment of the magnitude of the discretization error terms for the straight 
and curved element approximations can be obtained by considering a closed circular ring 
subjected to a harmonic lateral loading and approximated by either straight or curved 
elements. For this problem 

p = 0 (26) 

q = q Q sin ££ (27) 

If the exact solutions to equations (22), (23) and (24), (25) are denoted as u ex ,w ex and 
u ex> w ex> respectively, these two solutions are related by 



where p is the radius of gyration of the arch. 

Because the determination of the magnitude of the principal error terms of order 
h^ is quite tedious, selected numerical calculations were first carried out for the ring 
problem to investigate the magnitude of the error term. Results were obtained for sev- 
eral values of m and ring geometries. The numerical results indicated, as expected, 
that the convergence rate is proportional to 1/N2. Typical results for an arch with 
area of 1 inch^ (6.45 cm2), moment of inertia of 10 inch^ (416 cm^), radius of 
100 inches (2.54 m), and m = 2, approximated by 24 straight elements per half wave, 
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were e- = -0.001955 and e— = -0.000892 (see table n) where u = u ex ^l + e-j 
and w = w ex (l + e-). 

Results for the same problem based on a curved element approximation (using the 
same number of elements) were €u = e w = -0.136862. Because the size of e u and e w 
are quite large for the number of elements used for this example, the principal error 
terms of order h 2 were determined for the curved element approximation. These 
error terms for u = u ex (l + e u ) and w = w ex (l + e w ) (see table II) are 

e u = e w ^ x (m > 1) (30) 

12N 2 (£] (m 2 - l) 

Equation (30) shows that the principal error due to the curved element approximation is 
inversely proportional to the ratio (p/R) 2 , which for a real arch is small but finite. 
Although this particular curved element approximation converges to the correct result, 
it has significant error for calculations based on a moderate number of elements and is 
an undesirable one. Thus a curved element may be developed which converges to the 
right equation, but the error may still be substantially greater than for the straight ele- 
ment. The reason for the undesirable character of the error term in equation (30) and 
the accuracy of the curved element used herein may be due to the fact that the usual dis- 
placements u and w do not admit unstressed rigid body motion of the arch element. 
Reference 11 has shown with numerical examples that omission of the rigid body modes 
in the element behavior appears to affect significantly the accuracy of curved structure 
approximations . 


Two-Dimensional Plane Stress Elements 

The linear elastic plane stress equations for equilibrium in the x- and y-directions 
formulated in terms of displacements are, respectively 


1 - p 

u ’xx + 2 u, yy 


+ v > xy + ~ = o 


B 


( 31 ) 


1 + jU. 
2 


u 


’xy 


+ 


1 - n 
2 


v >xx 


Py 

+ v, yy + b = 


(32) 


Here u and v are the displacements in the x- and y-directions, respectively, p x and 

Ft A 

p are the distributed forces, respectively, p is Poisson's ratio, B = — is the 

y i - p 2 

extensional stiffness, and t is the plate thickness. Subscripts following a comma denote 
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differentiation with respect to the indicated variable. A satisfactory finite element 
approximation should lead to equations which converge to equations (31) and (32) at a 
node as the element size vanishes. 

Rectangular elements .- The error analysis procedure is extended to a general 
rectangular plane stress element in appendix C and subsequently specialized for two 
models the stiffness properties of which are documented in reference 12. For the linear 
stress model the stresses cr x and oy in the x- and y-directions, respectively, are 
assumed to vary linearly while the shear stress r X y is constant (fig. 4) . For the linear 
edge displacement model, the displacements along an edge of the element are assumed to 
vary linearly (fig. 4) . The nodal variables used to define the stiffness matrices for these 
finite elements are 



(33) 


and a typical finite element equation contains contributions from all elements contiguous 
to the node. The pattern arrangement composed of equal elements is also shown in fig- 
ure 4. The distributed forces on the plane stress body were again concentrated in a 
simple fashion based on the value of the distributed force at each node location. The 
typical finite element equations were obtained and the error terms evaluated. An inves- 
tigation of the convergence of the finite element equilibrium equation for both models 
showed them to converge to the plane stress equations (31) and (32) with a principal error 
of order h 2 . (See table 1(b) .) 

The principal error was evaluated for the case where the two types of rectangular 
elements were used to approximate a square plate in plane stress subjected to a harmonic 
inplane loading, in the y-direction. The loading is 

P x = 0 (34) 


„ ^ oi„ ni7rx . n7ry 

P.. = P„ sin — — sin — — = L 
*y a a 


(35) 


and the plate is supported on the boundary such that the force resultant in the x-direction 


and the displacement v 
in the displacements u 


both vanish. For ju = 0.3 and m = n, the errors e 
and v, respectively (see table II), are as follows: 


u 


and e 


v> 


Linear stress model 


N 2 


(36) 


2.40 

n 2 


( 37 ) 
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Linear edge displacement model 


e U = 


1.15 

N 2 


(38) 


e v = 


1,97 

N 2 


(39) 


where N is the number of elements per half wave length. 

Triangular elements .- Results were also obtained for the use of the classical 
triangular plate element (ref. 13) to approximate plane stress problems. Arrangements 
or patterns A, B, and C were investigated for the convergence of three right triangular 
elements (see fig. 5 for patterns about point i,j). It was found that pattern A converges to 
the required plane stress equations (see eqs. (31) and (32)) and the principal error is of 
order h^. However the corresponding equations for pattern B are 

u >xx + 4"^ u jyy + ^ = 0 ( 4 °) 

1 j iiv >xx + v ryy + ^ =0 ( 41 ) 

and for pattern C are 


U ’XX + 2 u, yy + ^ + ^ v ’xy + B " 0 

(42) 

1 p 

(1 + M)u, xy + 2 M v Jxx + v, yy + ^=0 

(43) 


o 

The additional error terms for all patterns proportional to h are given in table 1(b). 

Some points can be noted by comparing the convergence characteristics when 
h — 0 of the pattern B and C equations with equations (31) and (32). First of all both 
patterns B and C lead to a discretization error of order 1. The fact that the pattern B 
equations do not contain cross -derivative terms suggests that convergence of shear 
behavior at the nodal point is poor. This poor convergence is not unexpected since there 
is no mechanism in the finite element equations for representing changes in shear at the 
nodal point because of the arrangement of the elements. However the pattern C equa- 
tions overprescribe the cross -derivative term by a factor of 2. Note that the difficulty 
arises from the element arrangement rather than the element properties since the ele- 
ment used here fully represents all states of plane stress. These results indicate that 
convergence difficulties may arise for some triangular elements as a result of poor ele- 
ment arrangement even though the element is well formulated. 
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Since the difficulty with pattern B is due to its inability to represent the cross 
derivative at the node, better convergence properties would be expected if additional 
degrees of freedom were used to characterize the right triangular element behavior. 

Such added degrees of freedom might be the deflections at the midpoint of the various 
edges or the derivatives of displacements at nodes. 

Fortunately, if patterns B and C are used in structural idealizations, they usually 
occur in pairs (see fig. 5) and the underprediction of the shear stiffness at one point is 
compensated to some extent by an overprediction of shear stiffness at a neighboring point. 
Nevertheless these results suggest that caution should be exercised to ensure that an 
excessive number of either patterns B or C does not occur in a structure when the results 
are strongly dependent on shear stiffness. A more consistent approach is to use pat- 
tern A since it converges to the appropriate plane stress equation. 

For completeness the convergence of a pattern composed of equilateral triangles 
in plane stress was also investigated. The typical pattern is shown in figure 6 and the 
resulting finite element equations were found to converge to the plane stress equations 
with a principal error of order h^. (See table 1(b).) 


Two-Dimensional Plate Bending Elements 


The error analysis procedure was developed for a general rectangular plate bending 
element in appendix D and subsequently specialized for three models. The models inves- 
tigated were those developed by Papenfuss (ref. 14), Melosh (ref. 15), and one developed 
independently by Adini and Clough (see ref. 1 or ref. 16) and Melosh (ref. 17). (These 
will be denoted respectively as the Papenfuss, Melosh, and ACM models; the stiffness 
matrices of all three models are tabulated conveniently in ref. 1.) The pattern arrange- 
ment is shown in figure 7. The nodal variables for these finite element models are 


r ^ 



(44) 


where w is the lateral displacement and 6 and cp are the rotations about the y- and 
x-axes, respectively. On the basis of the beam results a consistent set of three plate 
bending finite element equations should be expected to converge to 



(45) 

e = w, x 

(46) 

II 

* 

<< 

(47) 
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as the element size vanishes. Here equation (45) is the familiar plate equation and equa- 
tions (46) and (47) are constraint equations between the rotations 6 and <p and the 
derivative of w. All three models lead to equations which converge to constraint equa- 
tions of the form of equations (46) and (47) and the Melosh and ACM models also converge 
to the plate equation (45). The Papenfuss model however converges to the equation 


4 

V w + 


2a; 2 2 2 

— + — - — + — ]w 

35 35qi2 25/ 


’xxyy D 


(48) 


where a is the aspect ratio of the element and thus has a principal error of order 1. 
The principal errors for all other equations for the three models are proportional to h 2 
with the exception of the constraint equations for the ACM model which are proportional 
to h 4 . (See table 1(c).) 

It is well known from numerical calculations (ref. 1) that the Papenfuss model has 
some deficiencies and that the source of the discrepancies is the inability of the model 
to describe the twist behavior of a plate. This discrepancy term shows up as an incor- 
rect cross -derivative term in equation (48) when h vanishes. 

Square elements of the three models were used to approximate a simply supported 
square plate subjected to a harmonic loading 


q = q G 


sin 


m7TX 

a 


sin 


n7TX 


(49) 


The same procedure was used to approximate the lateral vibration characteristics of the 
plate. The error in deflection and the error in frequency resulting from the 
finite element approximation for m = n and ju. = 0.3 (see table II) are as follows: 


Model 

e l 

e 2 

Papenfuss 

-0.0463 - °- 2646 
N 2 

-0.0486 - 0,2909 
N 2 

Melosh 

0.708 

0.708 


N 2 

N 2 

ACM 

1.069 

1.069 


N 2 

N 2 


where N = i s the number of elements per Fourier half wave. Figure 8 shows a 

h 

sketch of the plate and a plot of the ratio of the finite element results to the exact result 
for the three elements as 1/N 2 vanishes. The results show that in the limit as the size 
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of the element vanishes, the results for this problem based on the Papenfuss model have 
an error of approximately 5 percent. Note also that the Papenfuss model converges from 
below and the other two models converge from above. 

Since the Papenfuss model equations do not converge to the plate equation, a solu- 
tion was obtained for the resulting Papenfuss psuedoplate equation (48) for a rectangular 
planform subjected to a uniform pressure q, having simple support boundary conditions 
and approximated by elements having the same aspect ratio as the plate. The solution 
was obtained by means of classical Fourier series expansions of the load q and deflec- 
tion shape as 

5 = l <*> 

m=l n=l 

w = l f w mn sin S2E sl „ SSL (51) 

m=l n=l 


where a and b are the lengths of the plate in the x- and y-directions, respectively, 
and m and n are restricted to the odd integers. 


The quantities w mn are related to the pressure q by 

16q 


w mn = 


tt 3 D: 


mn 


m4 m 2 n 2 /52 2a 2 2 \ nf 

a^ a 2 b 2 3 ® 35 o’ 2 / b^ 


(52) 


The following table compares the center deflection w c of the Papenfuss psuedoplate 

w P D 


with the exact plate results for 


- 4 
qa 


•x 10 3 : 


a/b 

Papenfuss 

psuedoplate 

Exact plate 

1 

l/ 2 

3.87 

9.64 

4.06 

10.13 


These results indicate that the Papenfuss psuedoplate has an error of approximately 
5 percent for both cases. Numerical results obtained in reference 1 for these cases 
appear to be converging toward these analytical results. 
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Remarks on Finite Element and Finite Difference Approximations 

Some general results can be deduced by comparing the convergence and principal 
error results of the various element approximations. One result deals with the require- 
ment of interelement compatibility. Consideration of the elements for both plane stress 
and bending gave examples where this requirement was neither necessary nor sufficient 
to insure convergence. The rectangular linear stress model is not interelement com- 
patible in displacements and the rectangular Melosh and ACM bending models are not 
compatible in slope; yet these three converge to the correct equations. On the other hand 
the right triangular plane stress patterns B and C are interelement compatible in dis- 
placements and the Papenfuss plate bending element is compatible in slope and displace- 
ments and yet these elements and arrangements do not converge to the correct equations. 

Some influence of unequal length segments is seen from the change in principal 
errors for the bar approximations. For bar elements of equal length the principal error 
is proportional to h^, and for bars of unequal length, proportional to h. From the 
asymmetric character of the Taylor series expansion about the reference point similar 
reduction of the order of error would be expected for other elements. This slower rate 
of convergence suggests that results may be less accurate when structures are approxi- 
mated by unequal elements than when approximated by equal length segments. 

Comparison of the errors for the plane stress elements with those of the bar ele- 
ments for the case of harmonic loading indicates that approximately nine bar elements 
and 15 square plane stress elements per half wave in one direction are required for a 
1 -percent error. Approximately two beam elements and nine square Melosh plate 
bending elements are required for a 1 -percent error. Thus more elements are required 
per wavelength in one direction for two-dimensional behavior than for one -dimensional 
behavior to obtain the same degree of accuracy. This fact is important because practical 
complex structures such as stiffened plates or shells are two dimensional and usually are 
approximated by various combinations of one- and two-dimensional elements. Since the 
elements have varying degrees of accuracy, results obtained for a structure approximated 
by a combination of the elements may be biased in some sense rather than having uniform 
inaccuracies. 

Results are given in appendix E which compares finite element approximations and 
central finite difference approximations to differential equations for the same problem. 
The results show that for bar, beam, and plane stress equations some correlation can be 
made between finite difference equations and equations resulting from finite element 
approximations. The results for the harmonic load problem for plane stress show that 
neither method is always more accurate and that the value of Poisson's ratio can change 
the relative accuracies of finite difference and finite element results. The results also 
show that the finite element approximations to beam problems are usually more accurate 
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than central finite difference approximations because of the higher rate of convergence of 
the finite element approximation. 


CONCLUDING REMARKS 

* Basic data are presented on the convergence and accuracy of finite element equa- 
tions resulting from patterns and elements in common use. The elements studied include 
bar elements, beam elements, plane stress elements of rectangular and triangular shape, 
plate bending elements of rectangular shape and straight and curved arch elements. The 
results indicate that many of the elements and patterns have good convergence proper- 
ties; that is, the resulting finite element equations at a node converge to the continuum 
equations at the node as the element size vanishes. The results also indicate that some 
elements and/or patterns have poor convergence properties. Right triangular plane 
stress patterns, for example, can have undesirable convergence properties. It is shown 
that the requirement for interelement compatibility is neither necessary nor sufficient to 
guarantee convergence of the finite element equations. Results for arch approximations 
show that both straight elements and curved elements converge to the behavior of an arch 
structure; however for the problem considered, the straight element approximation gives 
substantially better results than the curved element derived herein. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., December 16, 1969. 
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APPENDIX A 


ONE -DIMENSIONAL BENDING ELEMENT 


When bending behavior is introduced, finite element models and the discretization 
error analysis procedure become more complex than for extensional elements. To indi- 
cate the additional features brought on by bending, consider a simple prismatic beam sub- 
jected to distributed load q and approximated by an assemblage of beam bending finite 
elements of equal length h (fig. 2) . The finite element nodal variables for this prob- 
lem are 



(Al) 


where w^ and are the displacement and rotation, respectively, at a node, and two 
finite element equilibrium equations are obtained at each node. These two equations are 
expanded in a Taylor series about the ith point (when w and 6 are considered as 
independent variables) to give 


-12w" - lrw iv - jj_ w vi 


h° 

1680 


w vm + 


. + 120 r +2h 2 0'” +^0 V 
10 


+ 55o eVil * 


-3iL 

" El 


(A2) 


6 = w' 



h 4 h 6 

— — w v + — — 
120 5040 


w vii + 




,6 


2160 


e vi 


(A3) 


where I is the moment of inertia of the beam cross section. Since beam behavior for a 
continuum is defined by only one independent variable w, it is useful to eliminate 6 
insofar as possible from the finite element equation. This elimination is done by differ- 
entiating equation (A3) to obtain expressions for derivatives of 8 and sequentially back 
substituting these derivatives into both equations (A2) and (A3). 


Equations (A2) and (A3) finally can be put in the form 


w 


IV _ 


720 


w vm + 


jL 

El 


(A4) 


8 - w’ + 


180 


wv + 


= 0 


(A 5) 


Equations (A4) and (A5) show that in the limit as h vanishes the beam finite element 
equations converge to the familiar beam equation and the correct constraint equation 
between the two finite element variables 8 and w. 
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APPENDIX B 


DERIVATION OF STIFFNESS MATRIX OF A CURVED BEAM 


In this section the stiffness matrix of a curved prismatic beam bending and exten- 
sional element is derived. The method used is the Rayleigh-Ritz method wherein the 
strain energy for the element is expressed in the matrix form 



(Bl) 


where is the set of generalized displacements that characterize the element and 

[K] is the stiffness matrix. The energy for a curved beam of radius R and length h 
based on the theory of reference 8 is 

U = II 0 h i'A® T ^® dAdS (B2) 


where 




R 2 

-R 

-R 2 ? 


-R 

1 

R? 


-R 2 £ 

R£ 

R 2 ? 2 


(B3) 


(B4) 


and 


e = w* 



(B5) 


In the preceding equations u and w are, respectively, the tangential and radial dis- 
placement of the element, and s and £ are, respectively, the coordinates along the 
length and through the depth of the element. 

The displacements are assumed to vary over the element length as 

w = a Q + ajS + a 2 s 2 + a 3 s 2 (B6) 


u = bQ + bjS 


(B7) 
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APPENDIX B - Continued 


where and are constants. The nodal variables may then be related to the 
constants as 


where 




0 

0 

1 

0 

h 

1 


0 

0 

0 

0 

h 2 

2h 



1 0 

0 0 

k 0 

1 h 

0 0 

1 h 

R R_ 





- 


W 1 

h 

u 2 

w 2 


A 




(B8) 


(B9) 


(BIO) 



(Bll) 


and the subscript 1 in the {b^ vector indicates the end of the element where s = 0, 
and the subscript 2 in the vector indicates the end of the element where s = h. Dif- 
ferentiation and substitution of equations (B6) and (B7) into equation (B3) yields 


{»} = !>]<?} 


(B12) 
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APPENDIX B - Continued 


where 

0 

0 

0 

0 

0 

1 


D*J - 

1 

s 

s 2 

s 2 

0 

0 

(B13) 
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i / 
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1 

R_ 



Equation (B2) through the use of equations (B12) and (B8) may be rewritten as 


U - | J 0 h j A { 6 > T [H-5 T [ B ] T [ C ][ B ] [h- 1 ]®^ ds (B14) 


Since [h] and [b] are independent of the integration over the cross section, the 
energy may be expressed as 


u = | I 0 h @ T [ H_1 ] T [ B ] T [ T ] ® ds (B15) 


where 



R(1 + Z) 



-(1 + Z) 

-r 2 z 


-(1 + Z) 

5* 1 + z > 

RZ 


-R 2 Z 

RZ 

R 3 Z 


The following integrals were used to obtain equations (B16) 


(B16) 


C dA 
JR-5 

r C dA 

JR-C 

C C 2 dA 
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^(1 + z ) 

ZA S 
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(B17) 


where Z is the Winkler-Bach constant for curved beams. Thus a comparison of equa- 
tions (Bl) and (B15) reveals that 



(B18) 
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APPENDIX B — Continued 


The elements of the stiffness matrix are 

Rll k 12 


M = 


k 12 k 13 

k 14 

k 15 

k 22 k 23 

k 24 

k 25 

CO 

CO 
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k 35 

Symmetric 
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(Equation continued on next page) 
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APPENDIX B - Concluded 


where 
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In the error analysis study in the body of the paper, it is assumed that £ « R. Then 
Z may be approximated by 

Z = -L- (B22) 

AR 2 
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APPENDIX C 


RECTANGULAR PLANE STRESS ELEMENTS 

This appendix presents the development of accuracy relationships for plane stress 
problems approximated by finite elements. Consider an isotropic constant thickness 
body in plane stress subjected to distributed inplane loadings p x and Py and approxi- 
mated by rectangular finite elements (fig. 4). It is presumed that the only generalized 
coordinates used to describe the element behavior are the displacements u and v at 
the nodes. The rectangular plane stress finite elements are connected at the corners; 
however the displacements of two adjacent elements need not match along a common 
boundary; that is, the elements need not be interelement compatible. 

The behavior of a body in plane stress is governed by two equations of equilibrium. 
In typical finite element approximations, two equilibrium equations also occur at each 
node if the stiffness matrix of a typical rectangular finite element has the general form 



where u^ and v^ are the displacements of the ith node in the x- and y-directions, 
respectively, and F Xi and F yi are the corresponding element nodal forces (fig. 4) , 
and the quantities kij are the stiffness matrix coefficients which depend on the behavior 
assumed within the element. The nodal concentrated load is taken to be the pressure 
evaluated at that grid point multiplied by 1/4 of the total area of the four elements sur- 
rounding the grid point. Use of equation (Cl) for each of four adjacent identical elements 
leads to the following two finite element equations at a node i,j 
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APPENDIX C - Concluded 


4k U\j + 2k 13(“i+l,j + "1-1 j) + k 15("i + l, S+ l + “l-lj+l * VlJ-1 +u i + l,j-l) 


+ 2k i7(»i, i+ i + «t j-i) + k i6(vi, j+ i - n-i j+i + v i-i,j-i - Vi,j-i) - p x “ h2 - 0 


(C2) 


k 25(“i + l, i+ l - "i-lj+l + "l-lj-l ' "l + l,j-l) + 4k 22 v i,j + 2k 24(v i+1;j + V, j) 

+ k 26( v i+l,)+l + Vl,j + 1 + T l-1 J-l * v i+l,j-l) + 2k 28( v l,j + l + v i,J-l) - Py“ h2 = 0 < C3 > 

The variables u and v in equations (C2) and (C3) are expanded in a two-dimensional 
Taylor series about the point i,j to give 

4 ( k ll + k 13 + k 15 + k 17) u + 2h2 [] k 13 + k 15) u ’xx + ® 2 ( k 15 + k 17) u >yy + 2ak 16 v >xj] 

+ + k 15) u >xxxx + ® a k 15 u ’xxyy + a ( k 15 + k 17) u >yyyy + 4k 16^ a ' v >xxxy + a2v >xyyy)J 


+ . . . - P x ah = 0 
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’xxyy 
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k 26) v ’xx + a2 ( k 26 + k 28) v »yy + 2ak 25 u >xy] 

+ a 4 ^2g + k- 2 %)V ,yyyy + 4k25^Q!U, XXX y + 0( 2 U, x yyy)j 

(C5) 


where the subscripts following a comma denote differentiation with respect to the indicated 
variable. Equations (C4) and (C5), which represent equilibrium in x- and y-directions, 
respectively, can be used to investigate the convergence characteristics of any plane 
stress rectangular element with a stiffness matrix of the form of equation (Cl). 
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APPENDIX D 


RECTANGULAR PLATE BENDING ELEMENTS 

The development of finite element models for plate bending problems is less 
obvious than for beams, and several plate bending models, each derived from different 
assumptions, are in current use. Many of the models have three independent variables 
at each grid point; namely, the displacement w, the rotation 9 about the y-axis, and 
the rotation cp about the x-axis. (See fig. 7.) This section describes the error analysis 
procedure for a rectangular isotropic constant -thickness plate bending element the 
behavior of which is described by these three nodal variables. For a rectangular model 
of length h in the x-direction and width a?h in the y-direction, the stiffness matrix can 
be put in the form 
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APPENDIX D - Continued 

where the quantities are the stiffness matrix coefficients which depend on the 

behavior assumed within the element and where the directions of rotations and moments 
are defined in figure 7. Let a continuous plate subjected to a lateral load q be repre- 
sented by an assemblage of rectangular finite elements the stiffness properties of which 
are given in the form of equation (Dl). The three finite equations at a node i,j are 


+ 2k 14(*i + lj +w 1 -lj) + » 17 (w lij+1 ♦Wi jJ . 1 ) + ki (10 (* 1+1J+1 


+ W. 
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= qatr 


(D2) 


^2p hk 22 e i,j + 2hk 25( S i + l,j + + 2hk 28 (»i |J+ l + w ‘2,ll(«i+l,]+l + *1-1, j+1 

+ *1-1 J-l * *1+1 J-l) + 2k 24(w 1+ l,j - * 1-1 j) + k 2,lo( w l+l,j+l - »i-l J+1 - »i-l,)-l 
+ w l+l,J-l) + “ hk 2,12(*’l+l,J + l - *1-1 J+1 + ^i-1,3-1 - *1+1, j-l)] = 0 < D3 > 

J^iahkjs^j +2ohk 36 ^ 1+1J + n-l,]) + 2 “ Ut 39(» , ij+l + *i,j-l) + “ Ut 3,12(«’i+l,J+l 
+ *1-1 J+1 + *1-1 J-l + *1+1, J-l) + 2k 37(w ljj+1 - Wij.x) + k 3,10(»i+l J+1 + «i-l,j+l 

" w i-l,)-l ■ w i+l j-l) + hk 3,n(*i+i,i+i ■ * 1-1 j+i + * 1 - 1 , )-i ■ e i+i, i-i)j = 0 <D4> 

where D is the bending stiffness of the plate. 

The quantities w, 6, and cp in equations (D3) and (D4) are expanded in the Taylor 
series about the point i,j and through a process of differentiation and back substitution the 
quantities 6 and cp are expressed in terms of w. The result is 
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Equation (D2) is similarly expanded in Taylor series and equations (D5) and (D6) are used 
to eliminate 8 and cp. The result is 
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(Equation continued on next page) 
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(Equation continued on next page) 
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APPENDIX D — Concluded 




(D9) 


Equation (D8) in terms of general stiffness coefficients represents the vertical equilibrium 
equation of a plate approximated by finite elements in terms of the displacement w, and 
equations (D5) and (D6) are the constraints between 6, cp, and w. Equations (D5), (D6), 
and (D8) can be used to investigate the convergence of the finite element equations at a 
typical node point when the stiffness properties of the element are given in the form of 
equation (Dl). 
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APPENDIX E 


COMPARISON OF FINITE DIFFERENCE AND 
FINITE ELEMENT APPROXIMATIONS 

It is interesting to compare finite element approximations with the widely used 
central difference approximations to structural differential equations of equilibrium. 
While correlation is not obvious for all classes of problems some observations can be 
made. 


One -Dimensional Extensional Bar 

The governing equation for a one -dimensional bar subjected to a distributed inplane 
load p is 


U,xx + EA 


where the subscripts following the comma indicate differentiation with respect to the 
indicated variable. Application of central finite difference approximations to equa- 
tion (El) for unequal nodal spacing h and ah gives 


(1 + a)h 


1 + a 


u i + a u i+l 


Comparison of equation (E2) with the finite element equation approximated by segments 
of unequal length (eq. (3)) indicates that for this simple problem the finite difference and 
finite element equations are identical approximations. 


Plane Stress 

The governing equilibrium equation in the x-direction for a body in plane stress is 
given by equation (31). The central finite difference approximation to equation (31) is 
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APPENDIX E - Continued 


The corresponding finite element equation using rectangular elements and the same 
nodal spacing is 
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where Cj takes the following values for the two finite models considered in this paper: 
Linear stress model 
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Linear edge displacement model 
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The Taylor series expansion of the finite element equations about the point i,j can be 
written as 
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Study of equation (E7) shows that the finite element and finite difference equations 
for equation (31) are related by 


E x (u,v) = D x (u,v) + Cjh 2 - - i — - 9 - 

1 Ax 2 Ay 2 


(E8) 


where E x (u,v) and D x (u,v) are the element and difference operators, respectively, 
and where the central finite difference operator for u >XX yy is 


A^u 


-2 




Ax 2 Ay 2 Q! 2 h^ 


<-* 


-2 > u 


-2 


(E9) 
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APPENDIX E - Continued 


A similar investigation of the y-direction equilibrium equation gives the relations 

between the element and difference equations for equation (32) 

\ 


Ey(u,v) = Dy(u,v) + c 



A 4 v 

Ax 2 Ay 2 


(E10) 


where Ey and Dy are the corresponding element and difference operators. Equa- 
tions (E8) and (E10) suggest that any set of consistent finite element or finite difference 
operators for a rectangular mesh where the variables are the displacements at the nodes 
are related by 


E x( u > v )i,j 


E y (u,v) 




= D. 




+ Ir 



= D y( u ’ v) i,j + h 


2 



A 4 u 5 A 4 v \ 

Ax 2 Ay 2 2 Ax 2 Ay 2 /^ 

A 4 u - A 4 v \ 

Ax 2 Ay 2 4 Ax 2 Ay 2 /^j 


(Ell) 


(E12) 


where the quantities c^ are constants. Usually as a result of symmetry in assumptions 
on behavior of u and v these quantities c^ will be related by 


c l = 


c 2 " 



(E13) 


If c i = 0, the element operator becomes the standard central finite difference operator. 

Since these approximations all have the same error of order h 2 , the element and 
difference equations converge at the same rate. Thus the relative accuracy of the various 
approximations will depend on the problem solved. If the twist terms u )XX yy and 
v >xxyy are ^ ar S e f° r a problem, the finite element answers will be less accurate than 
the finite difference answers. For example consider the problem treated previously by 
finite element approximation, namely, a square plane stress plate subjected to an inplane 
harmonic loading in the y-direction. A solution to this same problem by central finite 
differences gives the errors resulting from the finite difference approximations e u 
and e v as 


€ n — 


-h 2 ir 2 m 2 n 2 


24a 2 (1 - jn)(m2 


— 2y(ll - 2 m + 3/i 2 ) 


(El 4) 
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APPENDIX E - Continued 


, , - LlJL m 8 + 1 + 10 a_± H 2 m 6 n 2 4 3+ 23fJ_+ M 2 ,.-_3)l 3 . m 4 n 4 .-5 ± - V 2 - 3 M 3 m 2 n 6 il + 2 g_ ; „8 

_ ._hV__ZZ i § § 4 

6a^(l - m) m® + — ~ H m4n2 + (2 - n)m2n4 + 1 ~ H n® 


,(E15) 


The following table compares the error magnitude of the finite element results 
with finite difference results for m = n and for Poisson’s ratios // = 0 and 0.3: 



M 

= 0 

M = 

0.3 


e 

c u A 

e v N 2 

e u N 2 

e v N 2 

Finite difference 

-2.26 

0.21 

-3.13 

-0.67 

Linear stress 

2.06 

2.60 

1.85 

2.40 

Linear edge displacement 

1.44 

2.26 

1.15 

1.97 


The results indicate that neither the finite element nor the finite difference method is 
clearly superior for the approximate solution to plane stress problems. 


Beam Element 


When bending is introduced the correlation between finite difference and finite ele- 
ments is more complex. This complexity is caused by the introduction of rotation vari- 
ables at the nodes in the finite element procedure, whereas finite difference methods are 
not usually developed on that basis. Some correlation can however be made for the simple 
beam. Consider the beam equation 

w iv = 5- (E16) 


The fourth derivative at a point i can be approximated by using data only at the points 
i - 1, i, and i + 1 by passing a fifth-degree polynominal through these three points and by 
expressing this polynominal in terms of Wj and &i, where 6 j is the derivative of the 
function at the node point 


w(x) = w i + #iX + 



2w i + w i + l , g i-l - W \.,2 


4h 


~ 5w i _ 1 + 5w i+ i _ gj_i + 8g t + 9 1+1 \ _ a + / -w i _ 1 + 2w t - w i+1 
v 4h3 4h 2 / \ 2h 4 


0 : 


i-1 


e. 


4h 3 


i+1 'x 4 


+ 


3w j-l - 3w i+ l 
v 4h 5 


e i-l + 4e i + 0 i+l\ „5 
4h 4 / 


(E17) 
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APPENDIX E — Continued 


Differentiating w(x) four times and evaluating it at point i gives 

w i-l * 2w i - w l+l , ' 9 i-l * e i+l 
h 4 2h 3 

Substitution of equation (E18) into equation (E16) gives a finite difference equation which 
is identical with the finite element equation for vertical equilibrium. There remains 
however to relate the moment equilibrium finite element equation to finite differences. 
This is done by recognizing that the moment finite element equation is a constraint 
between slope and displacement variable. This side constraint between these variables 
can be written exactly as 

pi+1 

w i + i = w i-i + J l . 1 edx < E19) 



(E18) 


Approximating the integral in equation (E19) by Simpson's one -third rule gives 


w i + l = w i-l + 3pi-l +40 i + 0 i+l) 


(E20) 


Comparisons show that equation (E19) is identical with the moment equilibrium finite 
element equation. The implicit use of the Simpson's integration formula indicates why 
the finite element approximation has greater accuracy and a higher order of error, 
namely, because the error in Simpson's rule is of order h 3 , whereas central difference 
approximations are of order h 2 . Some comparison of the relative accuracies of the 
standard central difference and finite element approximations can be determined by con- 
sidering the simply supported beam example subjected to a harmonic loading. The cen- 
tral difference approximation to w* v is 

w i iV = £f( w i-2 - 4w i-l + 6w i - 4w i+l + w i+2) < E21 > 

which when expanded in a Taylor series and applied to the harmonic load example prob- 
lem is 

(E22) 

Solution to equation (E22) is 


w = w ex (l + 6 d) (E23) 

where w ex is the exact solution and the error due to the finite difference approxi- 
mations 


_ h 2 /iutA 2 _ 1 7T 2 

d 6 V L ) N d 2 6 


(E24) 
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APPENDIX E - Concluded 


where N<j is the number of finite difference segments. Note that the finite difference 
error is proportional to i/n^ whereas the finite element error is proportional to 
i/n 4 . (See eq. (19).) Thus the finite element approximation is more accurate than the 
finite difference approximation for the same element size. On the other hand, the ele- 
ment equations lead to twice as many variables as the difference equations. For a 
1 -percent error in this example approximately two finite element segments and 13 finite 
difference segments are required. 
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TABLE I.- FINITE ELEMENT DISCRETIZATION ERROR 







DO 


TABLE I.- FINITE ELEMENT DISCRETIZATION ERROR - Continued 


Variables 


Element 

Equation 

Rectangle, linear 

a) 

stress distribution 

(2) 

Rectangle, linear 

(1) 

edge displacement 

(2) 

Right triangle, 

(1) 

pattern A 

(2) 

Right triangle, 

U) 

pattern B 

(2) 

Right triangle, 

(1) 

pattern C 

(2) 


(1) 

Equilateral triangle 

(2) 


(b) Plane stress elements 

Governing differential equations 

(1) u >xx + i-=-ttu, yy + v, xy + 0 

( 2 ) — ^ u >xy + 1 2 ^ v, xx + v >yy + = 0 

Error terms (appearing in left-hand side of equations) 


+ l2\ U ’xxxx + 


+ 12 \ (1 + M)u,xxx y 


[|{1 - M) + a 2 ( 2 + M 2 )Ju )xxyy + “ 2u ’yyyy + (1 + ^ )v ’Xxxy + + ^ a 2 y >xyyyj> + 

xxxy + ^ + M-) a ^ u >X yyy + v >xxxx + f " v) + ( 2 + M 2 )j v >xxyy + a2v ’yyyyj 


+ J 2 j^’xxxx + “ /*) + 2a ^] u >xxyy + 1 \ ^ a2u >yyyy + ^ + ^) v >xxxy + U + ^ a2v »x: 

+ §{(1 + ^ u ’xxxy + U + ^^xyyy + v >xxxx + l? 2 ^ “ ^ + 2 ] v >xxyy + «% 

+ Y2^ U,XXXX + Qf 2 (H J:i ) u, yyyy + 2av, xxxy + 3a?2v >xxyy + 2ff3v ’xyyyJ + • • * 

+ T 2 ^ au,xxx y + 3a2u,xx yy + 2a3u >xyyy + ( 2 /i ) v, xxxx + a 2 v, yyyy] + • • * 

+ n( u ™ + H JL “ 2 u> yyyy) + * ■ • 


M)“ 2 Lmv + • • • 


- v,x y + 12 ( u,xxxx + a2u, yyyy) + * ■ * 
- u >xy + T2^i^ v,xxxx + ff2v, yyyy) + * * * 


+ nr iiv ’*y + i2 


jj^ixxxx + 


“ 2u >yyyy + 2 < 1 + 4 )v !xxxy + 2(1 + ii)ci 2 - 


L_ — ' ' 

. 1 t M u, xy + y^2(l + /x)u ?XXX y + 2(1 + p.) 0 * 2 u, xyyy + V JXXXX + a2v >yyyyj + • • • 


' v >xyyy 

' v »yyyy 


•y] + . . . 


1,2 ("1^ U,xxxx + u, xxyy + ^ 2 ^ u, yyyy + v, xxxy + v >xyyy) + 

h 2 (Hls^ u, xxxy + u, xyyy + 32 ^ v> xxxx + v >xxyy + “gg^ v >yyyy) + 




Element 


Melosh 


ACM 


Papenfuss 


TABLE I.- FINITE ELEMENT DISCRETIZATION ERROR - Concluded 
(c) Rectangular plate bending elements 

Governing differential equations 

(1) V 4 w = i 

(2) e = w, x 

(3) <p = w 


Equation 

(1) 

( 2 ) 

(3) 


(1) 

( 2 ) 

(3) 

( 1 ) 


( 2 ) 

(3) 



Error terms (appearing in left-hand side of equations) 

- JL + gffiw 

12 48 / ,xxxx yy 


0.2(1 . jl ) + al 

\6 12/ 48 


• 0(h2) 

■o(h2) 


+ h" 


* ») 


w ’xxxxyy + g ( 4 + 


>xxyyyy 


+ o(h 4 ) 

+ o(h 4 ) 


+ ^35 “ 2 + 350,2 + ^V’xxyy + h2 


( a 6 , 2 a ,4 274 a 2 61 146 

\3675 2625 91 875 2625 ~ 3675a 2 / xxxxyy 


( 1 ~ I- 2 274 

61 

146a 2 \ 

\3675a 6 2625a 4 91 875a 2 

“ 2625 ’ 

3675 ) ’xxyyyy 


-o(h2) 

■ o(h 2 ) 
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TABLE n.- ERROR TERMS FOR HARMONIC RESPONSE EXAMPLES 


Element 

Error terms 

Bar: 

_ __ mTTX 

P = P 0 sin — 

u = u ex (l + e) 

p o / L \ 2 . m 7 Tx 

\i 0 X — 1 ] sin — — — 

EA\m7r/ L 

“ 2 = "ex 2 ^ - £ ) 

f 

*2 

12N 2 ■ 

Beam: 

/v ~ m7TX 

q = q Q sin — 

w = w ex (l + e) 

w ex = — (-i-) 4 S in 
ex El \m ttJ l 

= ^ex 2 ^ ' £) 

“« 2 = 1 (tt ) 4 

e- ?f4 4 
720N 4 

\rch: 

P = o q = q 0 sin ^ 

Curved segments: 
u = u ex (l + £ u ) 

W = W ex (l +€ w ) 

-q R 4 

u ex = cos ^ 

f o \ 2 R 

EIm(m 2 - ij 

q ° R -in mx 

El(m 2 - i ) 2 R 

Curved segments: 

f Jr 2 

12N 2 ( m 2 - l) 2 (g 2 
€ - - * 2 (] 

12N 2 (m 2 - l ) 2 (^) 2 ' 


Straight segments: 


Straight segments: 




TABLE II.- ERROR TERMS FOR HARMONIC RESPONSE EXAMPLES - Concluded 


Element 


Error terms 


Plane stress rectangle: 
p =° 


p = p_ sin sin 
m 3 a a 

u = u ex (l + e„) 
v = v ex ( 1 +e T ) 

P„a 2 (l + p)mn 


Bjr 2 (l - p)(m 2 + n 2 ) 


cos 2^ cos 2JX 


2p a' 
o 


I ( m 2 + lzJ‘ n 2) 

A 2 L sin ”25 si, 


Bjt 2 (1 - p)( ro 2 + n2) 2 


siniSisin HSL 
a a 


Linear stress distribution: 

e _ h 2 n 2 m 2 n 2 (5 -2 p + p 2 ) 
“ ' 12a 2 (m 2 + n 2 ) 


2 2 lit m» + 13 - 16 ^ 3 ^ m 6 n 2 + “ - 26 ^ + V - 2 F 3 ± Jil m 4„4 + SjMj. tlU?Z ^l±jd m 2 „6 + (LULf n S 

hV_ J 4 4 4 \ 2 / 


6a2(l - p) 


Linear edge displacement: 

_ h 2 TT 2 m 2 n 2 (7 - lOp - p 2 ) 


n® + - - — m 4 n 2 + (2 - p)m 2 n^ + — ~ ^ n® 


24a 2 (m^ + n 2 )(l - p) 


hV 1 f JLm ' 

" “ m) 


8 + 11 - Ijtt - f* Z m 6 n 2 + m 4 n 4 + 11 - 25p -H3p 2 + p 2 m 2 n 6 + (L-jA%l 

2 . 8 Q i i L 


n® + ^ ~ P m^n 2 + (2 - p)m 2 n^ + - “ f 1 n® 


Rectangular plate: 

Melos h: 



q ~ q sin sin 

^ ^o a a 

12 a 2 ^ 

1.jl + mE 

,6 12 48, 

| m 2 n 2 
f m 2 + n 2 

w=w ex (i +€j) 

ACM: 






p(x) 



(b) Approximation of continuous bar by finite elements. 




F 2’ U 2 


(c) Forces and displacements in bar finite elements. 


^ v /I + Ct\ 

F i h 2 



B— f* 



(d) Equilibrium at ith grid point in finite element assemblage. 
Figure 1.- One-dimensional bar and its finite element representation. 
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Figure 2.- Beam element. 




( 1 = 0 - 

Straight element approximation 



Curved element approximation 


Figure 3.- Arch approximated by straight and curved elements. 
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Figure 6.- Equilateral triangular finite element pattern. 
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